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m ABSTRACT 

I A computational procedure suitable for the solution of equations of motion for multibody sys- 

tems is presented. The present procedure adopts a differential partitioning of the translational 
motions and the rotational motions. The translational equations of motion are then treated 
■ by either a conventional explicit or an implicit direct integration method. A principal feature of 

the present procedure is a nonli nearly implicit algorithm for updating rotations via the Euler 
^ four-parameter representation. The present procedure is applied to the rolling of a sphere 

B through a specified trajectory, which indicates that the procedure yields robust solutions. 

m 
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1. Introduction 


The numerical solution of the equations of motion for multibody systems has been contin- 
uously challenging the dynamist. In general, computer simulation of multibody dyn ami cal 
(MBD) systems requires a concerted integration of several computational aspects. These 
include a data structure for describing the system topology, computerized generation of 
the governing equations of motion, incorporation and accurate treatment of constraint 
conditions, implementation of suitable solution algorithms and easy interpretation of the 
simulation results. 

In the past, the task of formulating equations of motion has been a dominant concern to 
many dynamists (see, e.g., [1-6]). From the computational viewpoint, it can be said that 
the differences in existing formulations lie principally in their ways of incorporating con- 
straints [7-15] and in their resulting system topologies [2,16,17]. When the MBD systems 
become more complex, such as in the deployment of large space structures, streamlined 
accommodation of system topologies becomes a more important concern than elegance of 
formulation. This is because a flexible data structure can allow different modeling, different 
formulations and different solution techniques to be adapted to different subsystems. In a 
previous study [14-15], a new constraint treatment technique that can solve the constraint 
equations in a separate module from that for the translation and rotation variables was 
presented. A major feature of that study was to preserve the system topology for a variety 
of MBD systems. 

In order to provide a complete set of solution modules, the constraint solution module must 
be interfaced with a solution module for the primary variables, viz, the translational and 
rotational variables. It is generally agreed that the solution of the translational motions 
can be treated either by a conventional explicit or an implicit direct integration method. 

However, a wide spectrum of solution techniques has been proposed to integrate the ro- 
tational motions. Yet, along with constraint algorithms, many solution reliability and 
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efficiency issues in multibody simulation packages appear to hinge on how one solves the 
rotational degrees of freedom. This is especially true for flexible multibody systems wherein 
the higher frequency-response components are often due to rotational ocillatory motions. 
The objective of the present paper is thus to present a computational procedure for a 
robust and efficient treatment of rotational motions so that one can solve general MBD 
equations for a variety of system topologies. 

2. Equations of Motion for Multibody Systems 

The discrete equations of motion for flexible multibody systems can be expressed as [16]: 

Mu + D{u) + S(u) + BjfX N + B%\ H = f(t) (2 • 1) 

$jv(u,t)=0, $tf(u,t)=0 (2-2) 

where M is the mass matrix, D(-) is the generalized velocity-dependent force operator, S(-) 
is the internal force operator due to member flexibility, Bn and Bn are the gradients of the 
nonholonomic and holonomic constraints (2.2), A n and A# are the constraint forces, f(t) is 
the applied force, u is the generalized displacement vector, (') denotes time differentiation 
and ( ) T designates the matrix transposition. 

The solution of (2.1) and (2.2) consists of two tasks: the satisfaction of the constraint 
conditions (2.2) to obtain A and the computation of u from (2.1). Procedures to obtain 
A n and A# by satisfying (2.2) were presented in [14,15] and will be adopted in Numerical 
Experiments. Hence, we will concentrate on the computation of u. 

3. Differential Partitioning of the MBD Equations 

A basic difficulty in direct integration of (2.1) is that u and u are not directly integrable, 
except for some special kinematic configurations. This motivates us to partition u into the 
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translational acceleration vector, d, which is directly integrable and the angular accelera- 
tion vector, w, which is not, and treat them by a partitioned solution procedure [18-20], 
viz 


u = 




(3-1) 


The equations of motion (2.1) can be partitioned according to the above acceleration 
partitioning: 


where 


M d 0 
0 M. 





Dj(u) + S d (u) - BJA 1 
(li) -f Su,(u) - BlX j 


in which 


B = Bn -I- Bh, A = Xn + \h 


(3-2) 


(3-3) 


(3-4) 


To effect the node-by-node integration for the rotational degrees of freedom, we partition 
Ci further into 


U= [tt 1 , W 2 , J 


p I T 


(3-5) 


where is a (3x1) angular acceleration vector for the j-th node, 


= [u4 j) ,u4 j) ,u 4 j) j 


(3-6) 
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4. Review of Staggered Stabilized Procedure for Constraints 


For simplicity, we rewrite (2.1) and (2.2) as 

Mu + D(u) + S(u) + B t A = f 
$(u, u, t) = 0 

A penalty technique is introduced to obtain the constraints as follows 

A = j$(u, u, f), 0 < c << 1 

Differentiation of A with respect to time yields 

eA = flu + Bu + 

dt 

Substituting u from (4.1) into (4.4), we obtain 

eA + BM~ l B T X = BM~ x {t — D{ u) - S(u)) + Bu + ^ = r A 

dt 

Hence, the solution of (4.1) and (4.2) has been replaced by (3.2) and (4.5). 


(4-1) 

(4-2) 


(4-3) 


(4-4) 


(4-5) 


5. Nonlinearly Implicit Procedure for Large Rotations 
5.1 Euler Parameters and Angular Accelerations 

The four-parameter Euler representation of the angular velocity for each node is expressed 
as (see, e.g., Wittenburg[2l]): 


where 


q 


1 o 

2 ui 



q = A(w)q, 


q = L 9o <h <12 93 J 


T 


0 

—U>3 

U) 2 

U>3 

0 

-W 1 


w 1 

0 


<j) — [ U>2 ^3 J 


(5-1) 


(5-2) 
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and the nodal-designating superscript is omitted for notational simplicity. 


Time differentiation of (5.1) once more yields 

q = A(w) -q + A(w) -q (5-3) 

where A(co) is obtained by substituting Co for w in A(w). 

Note that (5.3) contains the constraint condition 

q T q = 1 (5-4) 

in its second-derivative form: 

q r q + q T q = 0 (5*5) 

5.2 Mid-Point Integration of Euler Parameters 

Suppose that we know the state variables, Co k and u k , at the k-th time step and we want to 
solve for q* +1 . Because of the special properties of A(w) and A(cj), one can take advantage 
of the following set of mid-point rules: 

r q*+* = q* + 6q k+ * 

< q fc+ $ = q fc -f 6q* + ^ (5 • 6) 

q fc+i = 2q fc+* - q fc, 6 = § 

where h is the stepsize. 

Substituting (5.1) and (5.3) into (5.6), one obtains 

f [I - $A(co fc+ *) - 5 2 A(u/* + *)]q fc+ * = [/ - 6A(<o k+ *) + £A(u; fc ))q fc , , ?) 

\[I -6A{u k +i)]q k+ * =q k +6A(<o k+ *)q k +*, (q*+*)^*+* = 0 

Since A{u k+ ^) and A(w fe+ ^) are not available, we approximate them by 

A(u; fc+ *) a A(«*) , A(w fc+ ^) cz A(u k ) (5-8) 
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which are equivalent to a Newton-like tangent approximation. Hence, we have from (5.7a) 
and (5.8) 

[I — 6A{(jj k ) — <5 2 A(w fc )] • q fc+ * = q fc (5-9) 

Once q fc+ * is obtained from (5.9), we can obtain q fc+ * from (5.7b) by 

[/ - 5A(w fc+ *)]q fc+ * = q fc + 6A(w* + *)q* + * (5 • 10) 


Finally, one can update w from the following formula: 

« = 2(?og - gg - 9og), g=|.9i 92 93 J, 9o+8 T g =1 (5-11) 

where g has the same form as d; as given by (5.2). Hence, if necessary, one can obtain u> 
from the equations of motion and iterate on q and q by (5.7). 


The skew-symmetric solution matrices in the above difference equations can be explicitly 


inverted via the formula: 


‘1 

—a 

-6 

— c 

-1 

' 1 

a 

b 

c 

a 

1 

c 

-b 

1 

—a 

1 

—c 

b 

b 

— c 

1 

a 

1 + a 2 + b 2 + c 2 

—b 

c 

1 

—a 

_ c 

6 

—a 

1 . 


—c 

—b 

a 

1 _ 


Hence, the solution of (5.9) and (5.10) becomes straightforward. 


(5 • 12) 


5.3 Update of New Angular Orientation 

Once q fc+l is computed from (5.9) and (5.6c), it is often required to compute the 
body-fixed basis vector, b = [hi b 2 baj in terms of the inertial basis vectors, 
e = [ ei ej e$ J . These two vectors are related by 

b = Re (5 • 13) 


where 


R = 


“2(9o + 9?) “I 

2(?i92 ~ 9o93) 

.2(9193 + 9092) 


2(9i 92 + 9093) 
2(9o + 92 ) “ 1 
2(9293 ~ 9o9i) 


2(9i93 - 9092)* 
2(9293 + 9o9i) 
2 (9o + 9s) - 1 . 


(5 • 14) 
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In order to satisfy the orthonormality of R, it is crucial to satisfy the two constraints, 
viz, q T q = 1 in computing q and q T q = 0 in computing q, respectively. This can be 
accomplished for q fc+ i by augmenting the constraint and solving the following equation 


by a Newton-like procedure: 


E, q 
q T 0 


Gl-ft} 


(5 • 15) 


where E, is the (4 x 4) solution matrix in the lefthand side of (5.7a) and b 9 is the righthand 
side vector of (5.7a), respectively. 


Similarly, the constraint q T q = 0 can be satisfied by 


q 

q T 0 




(5 • 16) 


where E, is the (4 x 4) solution matrix in the lefthand side of (5.10) and b q is the righthand 
side vector of (5.10), respectively. 

In addition, after the solution has converged at the (fc -I- |)-timestep, we must enforce 
the same two constraint conditions in updating q fc+l and q fc+1 . This is effected by the 
following simple procedures. 

To maintain the constraint q T q = 1 at the new (k+l)-th timestep, we employ 




(5 • 17) 


where I is (4 x 4) identity matrix. Similarly, to enforce q T q = 0, we use 




(5 • 18) 


This completes the present nonlinearly implicit procedure for large rotational motions. 
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6. Solution Procedures for MBD Equations 

We recall (3.2) and (4.5) to summarize the present overall computational procedures 


Mid = f d - Z> d (u) - 5 d (u) - f? d A (6-1) 

M„u, = f W -D„ (u) - 5 W (u) -BlX (6-2) 

eA + BM~ l B T X = r* = - D(u) -5(u)) + Bu+ ^ (6-3) 


The computational sequences are as follows. 


<r +i =cP-i+M" 

d" +1 = d" +M n+i 


( 6 - 4 ) 


' A(w*+ j) - A(w*) A(w* + *) =* A(w*) 

[/ - 5A(w*+*) - 5 2 A(w fc+ *)]q*+* = [ I - $A(w fc+ *) + 5A(w fc )]q fc 
< [I - 5A(u;*+i)]q* + * = q fc + iA(w*+*)q*+*, (q fc +*) r q fc+ * = 0 
q* +1 = 2q fc+ i - q* 

q* +1 = 2q fe+ * - q*, (q fe + 1 ) T q fc+1 = 0 
. w = 2(g 0 g - gg - ?og) 


It should be mentioned that the above procedures require only (4 x 4)-matrices which can 
be inverted explicitly via (5.12). 


(el + ^BM" 1 B r )A fc+ ^ = 5r* + * + eA* (6 • 6) 

We will now apply (6.4) - (6.6) together with (5.15) - (5.18) to some sample problems. 

7 . An Example — Dynamics of a Bowling Ball 

This problem was investigated by Huston et al [22], whose equations do not involve the 
constraint force, A. In the present analysis, we employ a formulation that incorporates the 
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constraint force as part of the system variables. Fig. 1 illustrates the ball with its radius 
a and an offset center r 0 that is to follow a sine curve 


y = stnx 


(7-1) 


The various matrices and vector quantities for (6.1)- (6.3) can be derived as 

*m 0 -mr 0 ei-b2 mroex-bx O' 

0 m -mroe 2 • b 2 mr 0 e 2 • bx 0 
M = Ji 0 0 

sym. J 2 0 

J 3 . 


1 0 — abj • e 2 — ab2 * 62 — ob 3 • e2 

B = 0 1 abi • ex ab2 • ei ab 3 • ei 

cosx —10 0 0 

D = -mr J WlW3ei ■ bl + w 2W3ex • b 2 - ( u> l + w|)ei • b 3 ) 
0 \ wiw 3 e2 • bx + W2W3e2 • b2 — (u/f + w|)e2 • b 3 J 

{ W2W3(*^2 — ^ 3 ) ^ 

Uf3Wl(J3-Jl) / 

WxU>2(«^X — *^ 2 ) J 

( e 3 b 2 ' 

f d = 0, f w = mgr 0 ^ -e 3 • bx ► 

l 0 ) 

d = L^jyJ 2, ,c0 = [wx u>2 w 3 J r ,A=[Ax A 2 A 3 J 

5(u) = 0, f = 0 


(7-2) 


(7-3) 


(7-4) 


(7-5) 


(7-6) 

(7-7) 


b = Re 


(7-8) 


w = 2(9og-gg-9og), g=[?x 92 93 J' 


(7-9) 


There is a total of eight variables in the foregoing equations of motion as given by (7.6). 
However, in adopting the present solution procedure — viz, (6.4)-(6.6) — we solve for nine 
variables. 

Numerical solutions of the rolling of sphere on a flat sinusoidal curve have been obtained 


with the data summarized in Table 1. 
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Table 1 Physical Dimensions and Initial Conditions for a Rolling Sphere 


m = 71.32JV, a = 10.9cm, tq = 0 or 0.15cm 
Ji = J 2 = J 3 = 2/5 ma 2 , e = 10 -6 
x° = y° = 0, {u >2 = ~u>i = 1, 0 J 3 = 0} 

x° = t/° = au>?, {?o = Mi = ?2 = ?3 = °> 

The ball track curve on the flat surface with time is shown in Fig. 2 for the case of no 
offset (ro = 0) and the corresponding angular velocities in Fig. 3. The time histories of 
the three constraint forces are shown in Fig. 4 wherein (Ai, A2) correspond to the x and 
y- component of the constraint force to maintain the rolling contact condition, and A3 is 
to maintain the sinusoidal trajectory as imposed by (7.1). Hence, the first two constraints 
are indicative of skidding phenomenon and the third corresponds to the steering force 
required in the ball manuevering. Notice that they exhibit highly nonlinear behavior while 
still periodic. 

The ball track projected on the ball itself is shown in Figs. 5 through 8, the ( x — z)-plane 
view in Fig. 5, the (y — z)-plane view in Fig. 6, the (z — y)-plane view in Fig. 7 and a 
three-dimensional trajectory in Fig. 8. 

Convergence studies have been performed with increasing stepsizes and it has been found 
that the present procedure, viz, (6.4)-(6.6), maintains both the solution accuracy and 
stability for the stepsize up to h < 0.15. Although not described herein, an implicit 
algorithm was tried out to solve the translational motions, d, thus replacing (6.4) by a 
corresponding implicit procedure. It was determined that, while larger stepsizes were 
permitted, we not only had to iterate to convergence at each time step but also required 
more computations due to matrix solutions at each step and/or iteration. We hope that 
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further studies will illuminate how one can profit by a combined use of explicit and implicit 
algorithms to solve the translational part of the equations of motion. 

Figures 9-11 show the y-direction surface traction of the ball with the offset center 
(ro = 0.15a), the angular velocities and the contraint forces. Note that, for the ball with 
no offset center, the average velocity in the y-direction is found from Fig. 2 to be about 
1.2 unit/sec. On the other hand, the corresponding average velocity with the offset center 
is about 0.35 unit /sec. from Fig. 9 or about 1/3 of the no-offset case. A more dramatic 
variation with the offset center ball is illustrated in its angular velocities as shown in Fig. 
10. Note that the angular velocities no longer exhibit periodic response, whereas they 
are periodic for the no-offset case ( see Fig. 3). Likewise, the steering force to follow the 
sinusoidal curve (y = sinx) becomes highly nonlinear, although nonlinearly periodic. The 
x and y-direction contact force to maintain the rolling contact condition between the ball 
and the surface, although bounded, manifests extremely nonlinear behavior. 

Our experience with the example problem indicates that the present computational pro- 
cedure for handling large rotational motions coupled with translational motions is robust 
and efficient. It is important to note that the present procedure traces not only the an- 
gular motions accurately but more importantly the constraint forces and the four Euler 
parameters (although these are not presented here). We hope to test the present com- 
putational procedure for larger and flexible structural systems in the coming months and 
report further results. 

8. Concluding Remarks 

A coputational procedure for an accurate and efficient solution of large singular motions has 
been presented. The present procedure treats each nodal angular orientation separately 
in terms of the four-parameter Euler formula. Hence, one deals with only (4 x 4)-matrices 
whose inversions are done explicitly. It has been found from our numerical experiments 
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that it is essential to enforce the constraint condition on the Euler parameters (5.4) and 
its time derivatives at each integration step. The present procedure is thus well suited to a 
partitioned solution procedure wherein the translational motions are integrated separately 
from that of the rotations motions. In addition, the present procedure interfaces easily 
with the constraint force solution package as discussed in detail in Park and Chiou[14,15], 

From a theoretical viewpoint, the present procedure can be considered a nonlinearly im- 
plicit procedure since the angular accelerations and angular velocities are used as agru- 
ments in the implicit solution matrices. The stability of the present procedure is at present 
difficult to assess and needs to be addressed in order for the procedure to take full advan- 
tage of its implicit characteristic. This stability question will require an involved analytical 
study and we hope to address this issue in the coming months. 

Acknowledgements 

The work reported herein was supported by NASA/Langley Research Center under Grant 
NAG - 1 - 756. The authors wish to thank Drs. Jerry Housner and Jeff Stroud for their 
keen interest and encouragement during the course of the present work. 


14 


References 


1. Hooker, W. and Margulies, G., “The Dynamical Attitude Equations for an N-body 
Satellite,” J. Astronautical Science , 12, 123-128 (1965). 

2. Roberson, R. and Wittenburg, J., “A Dynamical Formalism for an Arbitrary Number 
of Interconnected Rigid Bodies with Reference to the Problem of Satellite Attitude 
Control,” Proc. the Third Int. Congress of Automatic Control , Butterworth, London, 
1965. 

3. Likins, P., “Analytical Dynamics and Nonrigid Spacecraft Simulation,” Technical Re- 
port 32-1593, Jet Propulsion Laboratory, Pasadena, Ca., 1974. 

4. Ho, J., “The Direct Path Method for Deriving the Dynamic Equations of Motion of 
a Multibody Flexible Spacecraft with Topological Tree Configuraton,” AIAA Paper 
No. 74-786, AIAA Conference on Mechanics and Control of Flight, Anaheim, Calif., 
1974. 

5. De Veubeke, B. F., “The Dynamics of Flexible Bodies,” Int. J. Engng. Sci., 14, 1976, 
895-913. 

6. Kane, T. and Levinson, D., “Formulation of Equations of Motion for Complex Space- 
craft,” J. Guidance and Control , 3, 1980, 99-112. 

7. Walton, W. C. and Steeves, E. C., “A new matrix theorem and its application for es- 
tablishing independent coordinates for complex dynamical systems with constraints,” 
NASA TR-R326, 1969. 

8. Baumgarte, J. W., “Stabilization of constraints and integrals of motion in dynamical 
systems,” Comp. Meth. Appl. Mech. Engr., 1. 1-16 (1972). 

9. Baumgarte, J. W., “A New Method of Stabilization for Holonomic Constraints,” Jour- 
nal of Applied Mechanics , 50, 869-870 (1983). 

10. Lotstedt, P., “On a penalty function method for the simulation of mechanical sys- 
tems subject to constraints,” TRITA-NA-7919, 1979, Royal Institute of Technology, 
Stockholm, Sweden. 

11. Wehage, R. A. and Haug, E. J., “Generalized coordinate partitioning for dimension 
reduction in analysis of constrained dynamic systems,” ASME J. of Mech. Design, 
104, 1982, 247-255. 

12. Huston, R. L. and Kamman, J. W., “A discussion on constraint equations in multibody 
dynamics,” Mech. Res. Comm., 9, (1982) 251-256. 

13. Fuehrer, C. and Wallrapp, O., “A computer-oriented method for reducing linearized 
multibody system equations by incorporating constraints,” Comp. Meth. Appl. Mech. 
Engr., 46, (1984) 169-175. 


15 



14. Park, K. C., “Stabilization of Computational Procedures for Constrained Dynamical 
Systems: Formulation,” Presented at the 27th SDM Conf., San Antonio, Texas, May 
1986, Paper No. AIAA-86-0926-CP. 

15. Park, K. C. and Chiou, J. C., “Evaluation of Constraint Stabilization Procedures for 
Multibody Dynamical Systems,” Proc. the 28th Structures, Structural Dynamics and 
Materials Conf., Part 2A, 1987, Monterey, CA, AIAA Paper No. 87-0927, 769-773. 

16. Bodley, C. S., Devers, A. D., Park, A. C. and Frish, H. P., “A Digital Computer Pro- 
gram for the Dynamic Interaction Simulation of Controls and Structures (DISCOS),” 
NASA Technical Paper 1219, May 1978. 

17. Orlandea, N., Chase, M. A. and Calahan, D. A., “A sparsity-oriented approach to the 
dynamic analysis and design of mechanical systems - part I and II,” Trans. ASME, 
J. Eng. for Industry, Ser. B, 99, 1977, 773-784. 

18. Park, K. C., Felippa, C. A. and DeRuntz, J. A., “Stabilization of staggered solu- 
tion procedures for fluid-structure interaction analysis, ” in: Computational Methods 
for Fluid- Structure Interaction Problems, Belytschko, T. and Geers, T. L. (editors), 
ASME, AMD Vol. 26, New York, N. Y., 1977, 95-124. 

19. Park, K. C., “Partitioned Analysis Procedures for Coupled-Field Problems: Stability 
Analysis,” Journal of Applied Mechanics, 47, 1980, 370-378. 

20. Park, K. C. and Felippa, C. A., “Partitioned Analysis of Coupled Systems,” in: Com- 
putational Methods for Transient Analysis, T. Belytschko and T. J. R. Hughes (eds.), 
Elsevier Pub. Co., 1983, 157-219. 

21. Wittenburg, J., Dynamics of Systems of Rigid Bodies, B. G. Teubner, Stuttgart, 1977. 

22. Huston, R. L., Passerello, C., Winget, J. M. and Sears, J., “On the Dynamics of a 
Weighted Bowling Ball,” Journal of Applied Mechanics, 46, 1979, 937-943. 


16 



Ball Trajectory (y) 


2.000 

1.600 

1.200 

OAiA 

V ♦ uw 

0.400 

0.000 

- 0.400 

- 0.800 

- 1.200 

- 1.600 

- 2.000 





Angular Velocity (w) 
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Fig. 3 Angular velocities of the sphere with no offset 



Contstraint Forces(A) 
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Fig. 4 Histories of three constraint forces on the rolling sphere 
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Fig. 7 Ball track projected on the X - Y plane of sphere 







Fig. 8 Ball track projected on 3-D sphere surface 
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Fig. 9 Ball track on the flat surface with offset ( r 0 = 0.15a) 
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Coatstraint Forces(A) 
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Fig. 11 Time histories of constraint forces with offset center 



